# -------------------------------------------------------------------------------
"Functions to compute value functions and bond pricing kernel"
# -------------------------------------------------------------------------------

# -------------------------------------------------------------------------------
#####   SOLVER FOR STAGES 1 and 2 #####
# -------------------------------------------------------------------------------
function solve_main(ce::CE_Economy_M, EW_P::Array{Float64,3}, EW_I::Array{Float64,3}; penalty_Π=0.)
# -------------------------------------------------------------------------------

  #Create matrices to fill-in
    bpP_mat   = SharedArray{Float64,3}(ce.nb,ce.ny,ce.nz);
    bpI_mat   = SharedArray{Float64,3}(ce.nb,ce.ny,ce.nz);
    WR_P_mat  = SharedArray{Float64,3}(ce.nb,ce.ny,ce.nz);
    WR_I_mat  = SharedArray{Float64,3}(ce.nb,ce.ny,ce.nz);
    πI_mat    = SharedArray{Float64,3}(ce.nb,ce.ny,ce.nz);

  #Create interpolation objects for policies
    knots     = (ce.b_grid, ce.y_grid, ce.z_grid);
    Π_exp_int = interpolate(knots, ce.Π_exp,   Gridded(Linear()));
  #Compute and interpolate expected value (RHS of Bellman) [note: EW_j does not include transition]
    itp_EW_P  = interpolate(knots, ce.ΠT[1,1]*EW_P + ce.ΠT[1,2]*EW_I  ,     Gridded(Linear()));
    itp_EW_I  = interpolate(knots, ce.ΠT[2,2]*EW_I + ce.ΠT[2,1]*EW_P  ,     Gridded(Linear()));
  #Interpolate prices
    knots_f   = (ce.b_grid, ce.y_grid, ce.z_grid);
    itp_q     = interpolate(knots_f, ce.q_mat,   Gridded(Linear()));

@sync @parallel for i_z=1:ce.nz
z_tilde = ce.z_grid[i_z]
for (i_y, y_v)        in enumerate(ce.y_grid)
for (i_b, b_v)        in enumerate(ce.b_grid)

          # Conjectured misreport
          Π_exp = ce.Π_exp[i_b,i_y,i_z]

          #Posterior zeta_prime after message m [include transition]
          #-----------------------------------------------
          zeta_hat_L     =  F_fun(ce, Π_exp, 1,  z_tilde);        #m=1(=L)
          zeta_hat_NL    =  F_fun(ce, Π_exp, 2,  z_tilde);        #m=2(=NL)
          zeta_hat       =  [zeta_hat_L; zeta_hat_NL]
          z_prime_L      =  Markov_Switch_fx(ce, zeta_hat_L)
          z_prime_NL     =  Markov_Switch_fx(ce, zeta_hat_NL)
          z_prime        =  [z_prime_L; z_prime_NL]

          # ----------------------------------------------------------------------------
          # SOLVE FOR OPTIMAL BOND POLICIES:
          # ----------------------------------------------------------------------------
          fun_bp_fx(bp_v::Float64)  =          z_tilde    *     (    f(ce, 0.)       * fun_bp(ce, bp_v, 0.,     b_v, y_v, ce.Bval,  itp_q[bp_v, y_v, z_prime[1] ], itp_EW_P[bp_v, y_v, z_prime[1] ])     +
                                                                  (1-f(ce, 0.))      * fun_bp(ce, bp_v, 0.,     b_v, y_v, ce.Bval,  itp_q[bp_v, y_v, z_prime[2] ], itp_EW_P[bp_v, y_v, z_prime[2] ])  )  +
                                             (1-z_tilde)  *     (    f(ce, Π_exp)    * fun_bp(ce, bp_v, Π_exp,  b_v, y_v, ce.Bval,  itp_q[bp_v, y_v, z_prime[1] ], itp_EW_I[bp_v, y_v, z_prime[1] ])     +
                                                                  (1-f(ce, Π_exp))   * fun_bp(ce, bp_v, Π_exp,  b_v, y_v, ce.Bval,  itp_q[bp_v, y_v, z_prime[2] ], itp_EW_I[bp_v, y_v, z_prime[2] ])  );
          #Solver
          res_bp            =  optimize(fun_bp_fx, ce.b_grid[1], ce.b_grid[end], Brent(),  abs_tol=1e-5, iterations=300)
          bp_v                =  Optim.minimizer(res_bp)
          val_max             =  -Optim.minimum(res_bp)
            #Compare with lower bound solution
                tmp_val     = -fun_bp_fx(ce.b_grid[1])
                if tmp_val  > val_max
                val_max     = tmp_val
                bp_v        = ce.b_grid[1]
                end
            #Compare with previous solution
                bp_prev_sol = ce.bpP[i_b,i_y,i_z]
                tmp_val     = -fun_bp_fx(bp_prev_sol)
                if tmp_val  > val_max
                val_max     =  tmp_val
                bp_v        =  bp_prev_sol
                end
          #Compute VP matrix for P-type
          #--------------------------------------------------------
          WR_P_mat[ i_b, i_y, i_z]  =  -(    f(ce, 0.)      * fun_bp(ce, bp_v, 0.,     b_v, y_v, ce.Bval,  itp_q[bp_v, y_v, z_prime[1] ], itp_EW_P[bp_v, y_v, z_prime[1] ])    +
                                          (1-f(ce, 0.))     * fun_bp(ce, bp_v, 0.,     b_v, y_v, ce.Bval,  itp_q[bp_v, y_v, z_prime[2] ], itp_EW_P[bp_v, y_v, z_prime[2] ])    );

          #Store bond policies in main matrix (idem for both types)
          #--------------------------------------------------------
          bpP_mat[  i_b,i_y,i_z]     =  bp_v;
          bpI_mat[  i_b,i_y,i_z]     =  bp_v;

          # Optimal choice of π for the I-type and compute WR_I_mat
          #--------------------------------------------------------
          #Evaluate value functions at the new debt policy
          π_prev_sol                =  ce.Π_exp[i_b,i_y,i_z]
          qL                        =  itp_q[    bp_v, y_v, z_prime[1] ];
          qNL                       =  itp_q[    bp_v, y_v, z_prime[2] ];
          WI_exp_L                  =  itp_EW_I[bp_v, y_v, z_prime[1] ];
          WI_exp_NL                 =  itp_EW_I[bp_v, y_v, z_prime[2] ];
          fun_tmp_L(π_v::Float64)   =  fun_bp( ce, bp_v, π_v,    b_v, y_v, ce.Bval,  qL,     WI_exp_L );
          fun_tmp_NL(π_v::Float64)  =  fun_bp( ce, bp_v, π_v,    b_v, y_v, ce.Bval,  qNL,    WI_exp_NL);
          fun_tmp_π(π_v::Float64)   =  f(ce, π_v)*fun_tmp_L(π_v)  + (1.0-f(ce, π_v))*fun_tmp_NL(π_v)  + abs(π_v-π_prev_sol)*1e-5;
          resI             =  optimize(fun_tmp_π, ce.π_min, 0., Brent(), abs_tol=1e-5,  iterations=300)
          πpol             =  Optim.minimizer(resI)
          val_max          =  -Optim.minimum(resI)
          #Compare with previous solution [allow small tolerance]
             tmp_val       = -fun_tmp_π(π_prev_sol)
             if tmp_val    > (val_max-penalty_Π)
             val_max       = tmp_val
             πpol          = π_prev_sol
             end
          #Store solution for the I type
          πI_mat[   i_b, i_y, i_z]   =  πpol;
          WR_I_mat[ i_b, i_y, i_z]   =  val_max;
    end
    end
    end

  #STORE RESULTS:
    dist_bp_abs   = maximum(abs,ce.bpP-bpP_mat);
    dist_bp_norm  = norm(vec(ce.bpP-bpP_mat))^2;
    dist_π_norm   = norm(vec(ce.πpol-πI_mat))^2
    ce.bpP        = ce.damp * ce.bpP  + (1-ce.damp)*bpP_mat
    ce.bpI        = ce.damp * ce.bpI  + (1-ce.damp)*bpI_mat
    ce.πpol       = ce.damp * ce.πpol + (1-ce.damp)*πI_mat
    copy!(ce.WPr,  WR_P_mat);
    copy!(ce.WIr,  WR_I_mat);
return dist_bp_abs, dist_bp_norm, dist_π_norm
end
# -------------------------------------------------------------------------------


# -------------------------------------------------------------------------------
# Bond pricing kernel, q(.)
# -------------------------------------------------------------------------------
function compute_prices(ce::CE_Economy_M; IIB=0)

  # Create matrix to fill in
  q_new             = SharedArray{Float64,3}(ce.nb,ce.ny,ce.nz);
  # Interpolating objects
  knots             = (ce.b_grid, ce.y_grid, ce.z_grid)
  default_P_bel_int = interpolate(knots,  ce.default_P_bel,    Gridded(Linear()))   #States: (y,b, zeta_tilde)
  default_I_bel_int = interpolate(knots,  ce.default_I_bel,    Gridded(Linear()))
  Π_exp_int         = interpolate(knots,  ce.Π_exp,            Gridded(Linear()))
  bpP_exp_int       = interpolate(knots,  ce.bpP,              Gridded(Linear()))
  # Interpolate relevant prices
  knots_f           = (ce.b_grid, ce.y_grid, ce.z_grid)
  q_int             = interpolate(knots_f,  ce.q_mat,        Gridded(Linear()))   #States: (y,bp,zeta_prime)
  if IIB==1
  q_int             = interpolate(knots_f,  ce.q_IIB_mat,    Gridded(Linear()))
  end

  @sync @parallel for i_zp=1:ce.nz
  zp_v = ce.z_grid[i_zp]
  for (i_y, y_v)   in enumerate(ce.y_grid)
  for (i_bp, bp_v) in enumerate(ce.b_grid)

          #Initialize recursion
          q_Pv, q_Iv = 0.0, 0.0;

          #Next period output:
          for (i_yp, yp_v)   in enumerate(ce.y_grid)
          # Next-period default conjecture for P and I
          dfP_bel = default_P_bel_int[bp_v, yp_v, zp_v]
          dfI_bel = default_I_bel_int[bp_v, yp_v, zp_v]
          # Update of beliefs after observing d=0, given expected policies of default.
          zp_v_tilde  = F0_fun(ce,zp_v,0,dfI_bel,dfP_bel)
          # Conjectured future inflation misreport for I-type (subject to d'=0 is observed)
          Π_exp = Π_exp_int[bp_v,yp_v,zp_v_tilde]
          Π_exp = min(max(Π_exp,ce.π_min),0.)
          # Conjectured future bond policies (b'')
          bpp        = bpP_exp_int[bp_v, yp_v, zp_v_tilde]
          #Prob(m=L∣Π_exp)
          fpΠ   = f(ce,Π_exp)
          #Prob(m=L∣Π=0)=0
          fp0   = f(ce,0.)

          ### Next period message
          for mp=1:2  # 1=Lie; 2=NoLie
          #Probability of message m, given conjecture and type
          Prob_m_P   =  fp0    #m=L
          Prob_m_I   =  fpΠ    #m=L
          if mp==2
          Prob_m_P   =  1-fp0  #m=NL
          Prob_m_I   =  1-fpΠ  #m=NL
          end
          # Update Posterior if d=0 & m'=mp.
          ζp_hat    =  F_fun(ce, Π_exp, mp, zp_v_tilde)
          ζpp       =  Markov_Switch_fx(ce, ζp_hat )
          # Next-period price
          qp       =  q_int[bpp,yp_v,ζpp]

          #Update recursion
          if IIB==0
          q_Pv += ce.Πy[i_y, i_yp] * (1.0-dfP_bel) .*  Prob_m_P * ((1-ce.mb)*(ce.zb+qp )+ce.mb);
          q_Iv += ce.Πy[i_y, i_yp] * (1.0-dfI_bel) .*  Prob_m_I * ((1-ce.mb)*(ce.zb+qp )+ce.mb);
          else
          q_Pv += ce.Πy[i_y, i_yp] * (1.0-dfP_bel) .*  Prob_m_P * Repayment_fx_IIB(ce, 0.,    qp);
          q_Iv += ce.Πy[i_y, i_yp] * (1.0-dfI_bel) .*  Prob_m_I * Repayment_fx_IIB(ce, Π_exp, qp);
          end

         end #close mp
         end #close yp

  q_new[i_bp, i_y, i_zp] =  (zp_v*q_Pv + (1.0-zp_v)*q_Iv) / (1.0+ce.r)

  end
  end
  end

#Compute distances
dist_q, dist_q_abs = 0.,0.;

if IIB==0
  dist_q  	         = norm(vec( q_new - ce.q_mat)).^2;
  dist_q_abs         = maximum(abs, q_new - ce.q_mat);
  copy!(ce.q_mat, q_new)
else
  dist_q  	         = norm(vec( q_new - ce.q_IIB_mat)).^2;
  dist_q_abs         = maximum(abs, q_new - ce.q_IIB_mat);
  copy!(ce.q_IIB_mat, q_new)
end
return dist_q, dist_q_abs
end
# -----------------------------------------------------------------------------


####   STAGE 0 #####
# -------------------------------------------------------------------------
function compute_def_prob(ce::CE_Economy_M, Wr::Float64, Wd::Float64)
  # Assumes that Wd∼Wd*eps [i.e., small error to randomize default choice], with eps∼N(1,σ)
  μ             = 1.0
  σ             = ce.σ_d
  a             = Wr/Wd
  prob_def      = cdf(Normal(μ,σ), a)
  return prob_def
end
# --------------------------------------------------------------------------

# -------------------------------------------------------------------------------
function solve_default(ce::CE_Economy_M, EW_P::Array{Float64,3}, EW_I::Array{Float64,3}, EQ_P::Array{Float64,3}, EQ_I::Array{Float64,3})

#Create matrices to store results
  WP_new     = Array{Float64}(ce.nb, ce.ny, ce.nz)   #Value function
  WI_new     = Array{Float64}(ce.nb, ce.ny, ce.nz)
  WPd_new    = Array{Float64}(ce.nb, ce.ny, ce.nz)   #Value function if default today
  WId_new    = Array{Float64}(ce.nb, ce.ny, ce.nz)
  QP_new     = Array{Float64}(ce.nb, ce.ny, ce.nz)   #Value function if already in default
  QI_new     = Array{Float64}(ce.nb, ce.ny, ce.nz)
  default_I  = Array{Float64}(ce.nb, ce.ny, ce.nz)   #Optimal default choice
  default_P  = Array{Float64}(ce.nb, ce.ny, ce.nz)
#Define interpolation objects
  knots      = (ce.b_grid, ce.y_grid, ce.z_grid)
  WPr_int    = interpolate(knots, ce.WPr,  Gridded(Linear()));  # Value fx if repayment
  WIr_int    = interpolate(knots, ce.WIr,  Gridded(Linear()));
#Compute expected value (RHS of Bellman) [with Markov Transition]
  EW_P_int   = interpolate(knots, ce.ΠT[1,1]*EW_P + ce.ΠT[1,2]*EW_I  ,    Gridded(Linear()));   # Expected W if not in default
  EW_I_int   = interpolate(knots, ce.ΠT[2,2]*EW_I + ce.ΠT[2,1]*EW_P  ,    Gridded(Linear()));
  EQ_P_int   = interpolate(knots, ce.ΠT[1,1]*EQ_P + ce.ΠT[1,2]*EQ_I  ,    Gridded(Linear()));   # Expected W if already in default
  EQ_I_int   = interpolate(knots, ce.ΠT[2,2]*EQ_I + ce.ΠT[2,1]*EQ_P  ,    Gridded(Linear()));


for (i_z, z_v) in enumerate(ce.z_grid)
for (i_y, y_v) in enumerate(ce.y_grid)
for (i_b, b_v) in enumerate(ce.b_grid)

        # Posterior after observing d=0  [before inflation misreport & message m]
        ζ_tilde     =  F0_fun(ce, z_v, 0, ce.default_I_bel[i_b,i_y,i_z], ce.default_P_bel[i_b,i_y,i_z])
        # Posterior after observing d=1
        ζ_tilde_def =  F0_fun(ce,z_v, 1, ce.default_I_bel[i_b,i_y,i_z],  ce.default_P_bel[i_b,i_y,i_z])
        ζpD         =  Markov_Switch_fx(ce, ζ_tilde_def)

        # Value function of repayment
        WPr_val = WPr_int[b_v, y_v, ζ_tilde]
        WIr_val = WIr_int[b_v, y_v, ζ_tilde]

        # Value function of defaulting
        WPd_new[i_b, i_y, i_z] = ce.UdefP[i_y] + ce.β*( ce.θ*EW_P_int[0., y_v, ζpD] + (1.0-ce.θ)* EQ_P_int[0., y_v, ζpD] )
        WId_new[i_b, i_y, i_z] = ce.UdefI[i_y] + ce.β*( ce.θ*EW_I_int[0., y_v, ζpD] + (1.0-ce.θ)* EQ_I_int[0., y_v, ζpD] )

        # Value function if already in default
        ζpD                   = Markov_Switch_fx(ce, z_v)
        QP_new[i_b, i_y, i_z] = ce.UdefP[i_y] + ce.β*( ce.θ*EW_P_int[0., y_v, ζpD] + (1.0-ce.θ)* EQ_P_int[0., y_v, ζpD] )
        QI_new[i_b, i_y, i_z] = ce.UdefI[i_y] + ce.β*( ce.θ*EW_I_int[0., y_v, ζpD] + (1.0-ce.θ)* EQ_I_int[0., y_v, ζpD] )


  # Randomizing the default decision (small shock to value functions)
  # -------------------------------------------------------------------------
  def_prob_P               = compute_def_prob(ce, WPr_val, WPd_new[i_b, i_y, i_z] );
  def_prob_I               = compute_def_prob(ce, WIr_val, WId_new[i_b, i_y, i_z] );
  default_P[i_b, i_y, i_z] = def_prob_P
  default_I[i_b, i_y, i_z] = def_prob_I
  #Value Function (fx of y,b,z)
  WP_new[i_b, i_y, i_z]    = def_prob_P  * WPd_new[i_b, i_y, i_z]   +  (1-def_prob_P) * WPr_val;
  WI_new[i_b, i_y, i_z]    = def_prob_I  * WId_new[i_b, i_y, i_z]   +  (1-def_prob_I) * WIr_val;

end
end
end

#Compute Distances and update
# ------------------------------------------------------------------------
#Distance between value fx across iterations
  dist_W_norm  = max( norm(vec(    WI_new  - ce.WI)) ,  norm(vec(   WP_new  - ce.WP))   ).^2;
  dist_W_abs   = max( maximum(abs, WI_new  - ce.WI ),   maximum(abs, WP_new  - ce.WP )  );
#Distance between policy & beliefs
  dist_def_abs = max( maximum(abs,ce.default_P-default_P),  maximum(abs,ce.default_I-default_I)  );
#Store results
  copy!(ce.WP, WP_new)
  copy!(ce.WI, WI_new)
  copy!(ce.WPd, WPd_new)
  copy!(ce.WId, WId_new)
  copy!(ce.QP,QP_new)
  copy!(ce.QI,QI_new)
  copy!(ce.default_I, default_I)
  copy!(ce.default_P, default_P)
return dist_W_norm, dist_W_abs, dist_def_abs
end
# -------------------------------------------------------------------------


# -------------------------------------------------------------------------------
# Computing Expectations of v [does not include Markov switches]
# -------------------------------------------------------------------------------
function compute_EV(ce::CE_Economy_M, W::Array{Float64,3})
  EW = Array{Float64}(size(W))
  for (i_zp, zp_v)     in enumerate(ce.z_grid)
    for (i_y, y_v)     in enumerate(ce.y_grid)
      for (i_bp, bp_v) in enumerate(ce.b_grid)
          EV_tmp = 0.0
            for i_yp = 1:ce.ny
              EV_tmp += ce.Πy[i_y, i_yp]* W[i_bp, i_yp, i_zp]
            end
          EW[i_bp, i_y, i_zp] = EV_tmp
      end
    end
  end
  return EW
end
# -------------------------------------------------------------------------------



# --------------------------------------------------------------------------
function CONJECTURE_fx(ce::CE_Economy_M)
# -------------------------------------------------------------------------------
  #Default beliefs
  dist_def_bel_abs   = max(   maximum(abs,ce.default_I_bel-ce.default_I)    ,   maximum(abs,ce.default_P_bel-ce.default_P)   )
  dist_def_bel_norm  = max(   norm(vec((ce.default_I_bel-ce.default_I)))    ,   norm(vec((ce.default_P_bel-ce.default_P)))   ).^2
  #Misreport beliefs
  dist_Π_bel_abs     = maximum(abs, (ce.Π_exp - ce.πpol))
  dist_Π_bel_norm    = norm(vec((    ce.Π_exp - ce.πpol)))^2
  #Update default policies and misreport, with dampening. Store results
  tmp_I              =  ce.damp  * ce.default_I_bel + (1-ce.damp)  * ce.default_I
  tmp_P              =  ce.damp  * ce.default_P_bel + (1-ce.damp)  * ce.default_P
  tmp_π              =  ce.damp  * ce.Π_exp         + (1-ce.damp)  * ce.πpol
  copy!(ce.default_I_bel, tmp_I);
  copy!(ce.default_P_bel, tmp_P);
  copy!(ce.Π_exp, tmp_π);
return dist_def_bel_abs, dist_def_bel_norm,  dist_Π_bel_abs, dist_Π_bel_norm
end
# ---------------------------------------------------------------------------


# -------------------------------------------------------------------------------
# Prices in Secondary Markets
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
function compute_qA(ce::CE_Economy_M; IIB_val=0)

  #Interpolate objects
  knots      = (  ce.b_grid, ce.y_grid, ce.z_grid)
  Π_exp_int  = interpolate(knots,  ce.Π_exp, Gridded(Linear()))   #States: (y,b, zeta_tilde)
  bpP_int    = interpolate(knots,  ce.bpP,   Gridded(Linear()))   #States: (y,b, zeta_prime)
  if IIB_val==0
  q_mat      = ce.q_mat
  else
  q_mat      = ce.q_IIB_mat
  end
  knots_f    = (ce.b_grid, ce.y_grid, ce.z_grid)
  q_int      = interpolate(knots_f, q_mat,    Gridded(Linear()))   #States: (y,bp,zeta_prime)
  # Create matrices to fill in
  qA_new     = Array{Float64,3}(size(ce.q_mat))
  qB_new     = zeros(ce.nb, ce.ny,ce.nz,2)

  #Main loop across states
  for (i_z, z_tilde) in enumerate(ce.z_grid)
  for (i_y, y_v)     in enumerate(ce.y_grid)
  for (i_b, b_v)     in enumerate(ce.b_grid)

          #Initialize recursion
          qA_Pv, qA_Iv = 0.0, 0.0;
          #Conjectured inflation misreport for I-type (subject to d'=0 is observed)
          Π_exp = Π_exp_int[b_v, y_v, z_tilde]
          Π_exp = min(max(Π_exp,ce.π_min),0.)
          #Prob(m=L∣Π_exp)
          fpΠ   = f(ce,Π_exp)
          #Prob(m=L∣Π=0)=0
          fp0   = f(ce,0.)
          # Bond policy (b')
          bp    = bpP_int[b_v, y_v, z_tilde]

          # Next-period Message
          for mp=1:2                # 1=L; 2=NL
            #Probability of message m'
            Prob_m_P      =  fp0    #m=L
            Prob_m_I      =  fpΠ    #m=L
            if mp==2
            Prob_m_P      =  1-fp0  #m=NL
            Prob_m_I      =  1-fpΠ  #m=NL
            end
            #Updated Posterior if d'=0 & m'=mp.
            ζhat          = F_fun(ce, Π_exp, mp, z_tilde)
            ζp            = Markov_Switch_fx(ce,ζhat)
            #Next-period Bond Prices
            q_PM          = q_int[bp,y_v,ζp];
            qB,qB_P,qB_I  = qB_fx(ce, Π_exp, ζhat, q_PM; IIB=IIB_val);
            #Update the Recursion
            qA_Pv         += Prob_m_P * qB_P;
            qA_Iv         += Prob_m_I * qB_I;
            #Store results
            qB_new[i_b,i_y,i_z,mp] = qB
         end #close mp

   #Store results
   qA_new[i_b, i_y, i_z]    =  (z_tilde*qA_Pv + (1.0-z_tilde)*qA_Iv)

  end
  end
  end
return qA_new, qB_new
end
# -----------------------------------------------------------------------------

#------------------------------------------------------------------------------
function Iterate_qA(ce::CE_Economy_M)
#------------------------------------------------------------------------------
    #First iterate on q prices
    dist_q_abs, i=100.0, 1
      while i<100 && dist_q_abs>1e-6
      dist_q, dist_q_abs = compute_prices(ce, IIB=0);
      i=i+1
      end
    #Compute price for IIB [iteratre until convergence]
    copy!(ce.q_IIB_mat,ce.q_mat)
    dist_q_abs, i=100.0, 1
      while i<200 && dist_q_abs>1e-6
      dist_q, dist_q_abs = compute_prices(ce, IIB=1);
      i=i+1
      end
    #Compute qA for non-indexed debt
    #------------------------------------------------------------------------------
    ce.qA_mat, ce.qB_mat          = compute_qA(ce::CE_Economy_M, IIB_val=0)
    #Compute qA for IIB indexed debt
    #------------------------------------------------------------------------------
    ce.qA_IIB_mat, ce.qB_IIB_mat  = compute_qA(ce::CE_Economy_M, IIB_val=1)
Void
end
#--------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
